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1.0  INTRODUCTION 


We  address  here  methods  for  triangulating  orbital  frame  photography 
with  rigorous  accounting  for  spacecraft  dynamics.  Triangulation  of 
OH tuXaZ  photography  was  established  as  an  important  survey  tool,  by  the 
extensive  mapping  of  the  Moon  accomplished  by  triangulation  of  orbital 
photography  (beginning  with  the  Lunar  orbitor  photography  and  continuing 
through  the  Appolo  missions) , and  by  more  recent  extra-terrestial  and 
terrestial  applications  and  studies. 

Several  issues  need  to  be  focused  to  appreciate  the  developments 
herein.  Dynamical  models  of  the  camera  bearing  vehicle's  motion,  from 
one  view,  are  unnecessary.  As  has  been  the  tradition  in  triangulation 
of  aerial  photography  (and  mostly,  orbital  photography),  it  is  possible 
to  construct  workable  triangulation  algorithms  based  solely  upon  prin- 
ciples of  geometric  optics.  Perhaps  the  first  question  a skeptic  might 
raise  is:  "Why  clutter  the  (otherwise  algebraic  equation)  scene  with 
differential  equations  that  must  be  solved  numerically?"  Enforcement  of 
dynamical  constraints  imposes  a physical  truth  (to  an  essentially  negli- 
gible degree  of  approximation,  the  camera  motion  -C&  governed  exactly  by 
Newton's  laws  of  motion)  - as  a direct  consequence  of  introducing 
dynamical  constraints,  one  whould  expect  to  more  closely  recover  the 
true  camera  position.  The  number  of  unknown  parameters  subject  to 
differential  correction  is  usually  reduced  as  a consequence  of  intro- 
ducing dynamical  constraints;  thus  the  normal  equations  have  corres- 
pondingly reduced  dimensions.  Thus  computer  storage,  run  time,  and 
precision  considerations  generally  favor  the  dynamically  constrained 
approach.  The  use  of  a dynamical  model  in  orbital  triangulation  allows 
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one  to  restrucuture  the  triangulation  process  as  an  "iterative,  extented 
Kalman  estimation  algorithm"  with  most  significant  computational  ad- 


vantages. 


The  above  comments  are  fully  supported  by  the  results  presented 


herein. 


2.0  DYNAMICALLY  CONSTRAINED  TRIANGULATION:  SUMMARY  OF  BASIC  EQUATIONS 

2.1  Colinearity  Equations 

In  the  process  of  satellite  photogrammetric  triangulation,  the 
earth  fixed  object  space  coordinates  of  various  earth  surface  features 
can  be  deduced  from  satellite  photograph  coordinates  of  the  images  of 
those  features.  The  fundamental  mathematical  transformation  central  to 
this  process  is  based  solely  on  the  principles  of  geometric  optics  and 
can  be  deduced  from  the  geometry  shown  in  Figures  1 and  2.  This  trans- 
formation (the  colinearity  equations)  can  be  written  [1]  as 


x - f 
o 


- f 


pnW  + ci2<w  + ci3(zp-zc)"| 
LC31<VV  + C32<VV  + C33<VVJ 

rc21(y*c>  + C22(VV  + C23(VZc>~] 
LC31(VV  + C32(VV  + C33(VZc>J 


Equations  (1)  project  the  position  of  a point  located  at  objtcX  coordi- 
nates (Xp,  Yp,  Zp)  into  its  -image  plane  coordinates  (x^,  y ) for  a 
camera  with  principal  point  offset  (xq,  y^)  and  focal  length  f.  The 
perspective  center  of  the  camera  is  located  at  objzcX  space  coordinates 
(Xc,  Y^,  Z^) . The  AjnCLQQ.  space  (x,y,z)  axes  are  oriented  relative  to  the 
objzcX  space  (X,Y,Z)  axes  by  the  direction  cosines  * 1,2,3. 

Equations  (1)  can  be  written  functionally  as 
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(2) 


x = F (X  ,Y  ,Z  ,X  ,Y  ,Z  ,C.  ,x  ,f) 
p p*  p*  p‘  c’  c’  c’  *m.  o’ 

y = G(X  ,Y  ,Z  ,X  ,Y  ,Z  ,C„  ,y  ,f). 

P P P P C C X *mj  ° 

If  several  points  are  considered  in  each  of  n photographs,  Equations  (1) 
and  (2)  can  be  doubly  subscripted  to  denote  the  equations  corresponding 
to  the  ith  point  in  the  jth  photograph  as 


t = F (X  ,Y  ,Z  ,X  ,Y  ,Z  ,C.  , x ,f.) 

"u  ^ 'I  Pi  cj  V J 


{ = G(X  ,Y  ,Z  ,X  ,Y  ,Z  ,C„  ,y  ,fj), 

pij  Pi'  pi’  pi'  'i'  p/  'j'  ‘-j'S’ 


(3) 


Equations  (3)  can  be  segregated  into  subsets  depending  on  whether  or  not 
the  particular  point  is  a contAoZ  poZnt  (i.e.  one  for  which  there  exists 
apriori  knowledge  of  its  object  space  coordinates).  Other  points  having 
distinctly  measurable  image  plane  coordinates  in  two  or  more  photographs 
are  unknown  and  are  called  pa6A  point* . Depending  upon  the  number  of 
photographs  in  a strip,  the  number  of  strips  with  sidelap,  the  number 
and  distribution  of  control  points,  the  number  of  pass  points  and  the 
manner  in  which  they  are  shared  by  adjacent  photographs,  an  intimidating 
variety  of  observation  equations  of  the  form  (3)  can  be  defined  in 
practical  applications.  It  is  not  uncommon  to  encounter  several  thou- 
sand such  nonlinear  equations  containing  several  hundred  unknowns.  An 
obvious  computational  burden  is  associated  with  the  least  squares 
solution  of  this  class  of  problems. 

In  the  traditional  approach  to  this  problem,  the  camera  center 

coordinates  (X  ,Y  ,Z  ) and  three  camera  orientation  angles 

Cj  Cj  Cj  J J J 

for  each  photograph  (it  is  possible  to  replace  the  direction  cosines 

C^;  as  functions  of  three  orientation  angles),  as  well  as  the  pass 

point  coordinates  have  been  treated  as  independent  unknowns.  The 

introduction  of  orbital  dynamical  constraints,  conceived  by  Brown  [2,3] 
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the  elimination  of  the  3n  (n  = number  of  photographs)  camera  center 
coordinates  in  favor  of  6 osculating  orbital  elements.  One  purpose  of 
this  work  is  to  carry  the  incorporation  of  dynamical  constraints  to  its 
logical  conclusion  through  rigorous  satisfaction  of  the  satellite 
equations  of  rotational  motion.  This  process  allows  the  further  re- 
duction of  unknowns  from  3n  orientation  angles  to  6 osculating  attitude 
elements  or  constants  of  the  rotational  motion. 

2 . 2 Orbital  Dynamical  Constraints 

The  essence  of  orbit  constrained  photogrammetry  is  the  recognition 
that  the  camera  exposure  stations  along  a given  strip  of  n photographs 
are  dynamically  constrained  according  to 

Xc  = Xc^ti,cl,c2’"',c6^ 
j J 

Y — Y(t,,c_,c->...fc^.)  j — l»2>*..>n,  (^) 

c_.  c J i i.  o 

ZCj  = Zc(tj’c1*c2 c6) 

where  Equations  (4)  are  functional  representations  of  the  solution  of 
the  spacecraft's  translational  equations  of  motion.  The  6 constants 
(cis...,c.)  can  be  any  set  of  initial  conditions  or  osculating  orbital 
elements  which  uniquely  define  the  orbit.  The  least  squares  solution  of 
the  photogrammetry  problem  is  modified  to  include  the  dynamical  con- 
straints as  follows:  Current  estimates  of  c, ,c0,...,  c,  are  used  in 
Equations  (4)  to  compute  X ,Y  ,Z  at  each  photograph  exposure  time 

Cj  Cj  Cj 

tj , the  resulting  estimates  of  camera  exposure  coordinates  are  used  in 
Equations  (3)  together  with  the  other  parameter  estimates  (X^  ,Y^  ,Z^  , 
4>.  ,0. ,y  ,f  ) to  determine  current  computed  values  of  x ,y 

J J J °j  °j  J Pij  Pij 

which  are  in  turn  dirrerenced  from  the  corresponding  observed  values  to 
find  the  current  residual  vector.  The  partial  derivatives  of  the 
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observation  Equations  (3)  with  respect  to  the  orbital  elements,  c^. 


1 


c_,...,c,  (the  elements  of  the  "A"  matrix)  are  determined  by  the  chain 
Z O 

differentiation  rule  applied  to  Equations  (3)  and  (4).  For  example 

3x  3X  3Y  3Z 

Pij  = 3F  °j_  3F  3F  

3c,  3X  3c,  3Y  3c,  3Z  3c,  * 

1 c . 1 c . 1 c . 1 

1 1 J 

The  partial  derivatives  with  respect  to  the  other  parameters  are  iden- 
tical in  form  to  the  corresponding  equations  for  unconstrained  photo- 
grammetry.  Depending  on  problem  requirements,  various  specific  forms  of 
Equations  (4)  are  available  and  the  corresponding  partial  derivatives 
readily  obtainable.  Further  discussion  of  orbital  constraints  is  found 
in  the  work  of  Brown  [2,3],  Light  [4],  Hartwell  [5],  and  Blanton  and 
Junkins  [8]. 

As  an  alternative  to  recovery  of  the  orbit  simultaneously  with  the 
triangulation  process,  one  can  use  tracking  data  to  recover  the  orbit 
apAtoAi.,  to  within  small  clock  biases,  and  often  acceptable  precision. 
The  photomeasurements,  in  most  cases  will  provide  only  marginal  improve- 
ment in  the  orbit  resolution  due  to  the  recent  advances  made  in  tracking 
systems  for  orbit  determination.  In  any  event,  the  triangulation  should 
be  dynamically  constrained  in  the  sense  that  the  sequence  of  exposures 
lie  along  a dynamical  path. 

2.2  Rotational  Dynamical  Constraints 

The  incorporation  of  rotational  dynamical  constraints  is,  in  con- 
cept, very  analogous  to  introduction  orbital  constraints.  The  details, 
however,  are  somewhat  more  complicated.  The  direction  cosines 
contained  in  Equation  (1)  can  be  written  as  functions  of  three  angles 
orienting  the  camera  {cj  axes  relative  to  the  earth  {ej  axes  as 

{c}  - [C]{e}  (5) 

where 
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This  is  the  classical  approach;  each  set  of  three  orientation  angles  (for 
each  of  n exposures)  are  treated  as  n independent  sets  of  three  unknowns 
in  the  triangulation  process. 


Perhaps  the  most  obvious  approach  (to  introduce  rotational  dynamical 

constraints)  is  to  derive  a set  of  differential  equations  describing  the 

(^ce’^ce’^ce^  an8les  and  recover  the  six  initial  conditions  (on  these 

three  angles  and  their  first  time  derivatives)  in  lieu  to  treating  the 

3n  angles  as  independent  unknowns.  Several  factors  should  be  considered, 

I however.  First,  rotational  dynamics  is  most  naturally  approached  by 

viewing  the  motion  from  a non-rotating  reference  frame  rather  than  an 

earth-fixed  frame.  Second,  care  should  be  taken  to  avoid  the  singularity 

that  exists  for  certain  orientations  [0  = + tt/2,  in  Eqn.  (6)]  for  any 

choice  of  orientation  angles.  One  should  allow  efficient  advantage  to 

be  taken  of  on-board  angular  motion  measurements  (e.g.  rate  gyros,  star 

sensors,  sun  sensors,  etc.).  In  view  of  these  considerations,  we 

introduce  (in  lieu  of  4>  ,0  , ) the  angles  (4>  ,0  .'i'  ) orienting  the 

ce  ce  Tce  co  co  co 

camera  relative  to  an  "orbiting"  frame  which  maintains  axes  in  the  local 
radial,  transverse,  and  normal  directions  so  that  (see  Appendix  2) 

[C]  = [CN]  = [ CO ] [ 00 ’ ] [ 0 ’ N ] [ EN ] T (7) 

} I where 


For  rotational  compaction,  we  use 


n 


C0 

co 


c"  for 


"cosine"  and  "s"  for 


(8) 


"sine". 
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where 


/ 0)  \ 
i x) 

/ 

1 wy/ 

= [CG]  | 

[ 0)  1 
\ Z / 

1^  = spacecraft  moments  of  inertia  about  mass  center  with  axes 
oriented  parallel  to  image  plane  axes. 

fc)  ,u>  ,u>  } = spacecraft  angular  velocity,  components  along  camera 
X ^ Z image  plane  axes 

,L  } = effective  external  torque  acting  on  the  spacecraft; 

components  taken  about  mass  center  with  axes  parallel 
to  camera  image  plane  axes. 

to  i 
gx 

<0  ( 
gy 

0)  l 
gz } 

(w  ,10  ,uj  ) = angular  velocity  components  along  rata  gyro's  three 
gx  gy  gz  orthogonal  axes. 

[CG]  = assumed  constant  3x3  interlock  direction  cosine  matrix, 
orienting  the  camera  axes  with  respect  to  the  rate  gyro 
axes . 

In  the  event  that  the  particular  spacecraft  departs  significantly  from  the 
rigid  body  assumption  (e.g.,  significant  internal  moving  masses,  vibration 
of  flexible  appendages),  additional  forcing  terms  will  appear  on  the  RHS 
of  Eqn.  (13),  and  one  additional,  coupled,  differential  equation  for  each 
such  mechanical  degree  of  freedom  must  be  written.  The  kinematical 
Equation  (13),  however,  holds  rigorously  for  the  motion  of  the  camera  axes, 
regardless  of  what  torques  are  present  and  regardless  of  departures  of 
the  actual  body  model  from  the  rigid  body  assumption.  To  employ  Eqns . (12) 
as  the  starting  point  for  modeling  rotational  dynamics  requires  excellent 
mathematical  models  for  the  actual  torques  (L^.L^.L^)  present  and  for 
any  significant  departures  from  the  rigid  body  assumption.  This  viewpoint 
was  pursued  initially  under  this  contract  and  yielded  in  the  results 
published  in  References  7 and  8 (abstracted  below  in  section  3.1). 

The  fact  that  Eqn.  (13)  is  rigorously  valid  (regardless  of  the  torque 
history  and  actual  body  flexibility)  and  since  true  angular  velocity  can  be 
measured  (by  conventional  rate  gyros)  motivates  a (tiJiteA  HoXt  i.nA.e.g>Ultion 
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approach  to  imposing  the  rotational  motion  constraint.  This  approach 
leads  one  to  consider  the  wisdom  of  bypassing  the  necessity  of  modeling 
the  spacecraft  flexibility  and  torque  history  by  directly  integrating 
Equations  (13).  This  latter  viewpoint  has  been  investigated  during  this 
research  project  and  doe.;  prove  to  be  the  most  attractive  option. 

The  orthogonal  axes  of  the  rate  gyro  package  generally  has  some 
orientation  other  than  the  camera  axes,  thus  we  must  account  for  the 
gyro-to-camera  interlock  rotation  as 


(14) 


where  [CG]  is  a nominally  constant,  3x3  interlock  direction  cosine 
matrix,  and  the  measured  angular  rates  are 


/ u>  \ 

to  V 

gxl 

(Ui  \ = [CG] 

'to  l 

) y( 

i 

f 00  ) I 

to  ] 

\ z 1 

gz/ 

to 

gx 

\ 

(ill  1 

l gxl 

/ bx 

(d  ; 

gx 

i 

)8yl 

+ /by 

<l> 

gz 

) 

(ill  1 

l gzl 

Ibz 

(15) 


{Measured}  = {True}  + {Bias}  +{Measurement} 

Noise 

If  one  substitutes  the  raw  data  (w  ,to  ,w  ) through  (14)  directly  into 

gx’  gy  gz 

(13)  and  integrates  to  determine  the  angular  rate  history,  the  effects  of 

the  uncompensated  rate  biases  {b  ,b  ,b  } and  noise  (v  ,v  ,v  ) will 

x y z x y z 

typically  result  in  unacceptably  rapid  divergence  of  the  integrated  orien- 
tation history  from  the  true  orientation  history.  The  results  of  numerical 
experimentation  with  typical  noise  levels  (1  arc  sec/sec)  and  bias  values 
led  us  to  conclude  that  the  measured  rates  should  be  passed  through  a 
noise  pre-filter  (section  4)  to  obtain  a filtered  angular  velocity 
history: 
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rf-  [C]  = [■£-  [CO]]  [00’]  [O’N]  [EN]T 
dPi 


(18) 


where  we  need 


a 3*rn(t) 

9 [CO]  = co 


3p 


3P-, 


90  (t) 

co 


3P, 


dip  (t) 

CO 


3P, 


“C-* 

CO 

sip 

Tco 

o' 

c0 

CO 

0 

-S0 

co 

n 

0 

0 

-s  ip 

CO 

Clp 

CO 

0 

0 

1 

0 

0 

-sd> 

CO 

C<f> 

CO 

0 

0 

1 

C0 

1—  CO 

0 

C0 

CO-J 

0 

-C<J> 

CO 

-s<i 

Tco  — 1 

CO 

sip 

rco 

o- 

-s0 

CO 

0 

-C0 

CO 

I 

0 

0 

-sip 

CO 

Clp 

TCO 

0 

0 

0 

0 

0 

c4> 

CO 

s<t> 

CO 

0 

0 

1_ 

C0 

— CO 

0 

-S0 

co  — 

0 

-sd) 

CO 

cd) 

CO  — 1 





— 

-sip 

CO 

cip 

CO 

0 

C0 

CO 

0 

-S0 

CO 

1 

0 

0 

— Clp 

CO 

— sUi 

co 

0 

0 

l 

0 

0 

c.4> 

CO 

sd> 

CO 

_ 0 

0 

0_ 

S0 

L-  CO 

0 

C0 

CO  ” 

0 

— S<J> 

CO 

c d>  _ 

CO 

1,2,3. .6,  (19) 


The  partials  [ 


3[^o»(t)’9co(t)’*co(t:)] 

3P, 


] are  obtained  by  integrating  the 


matrix  differential  equations 


$(t,tQ)  = F(t)4>(t,tQ) 


10  0 
0 10 
0 0 1 


and 


4-(t,to)  = F(t)4'(t,tQ)  + B(t) 


f(to’to)  = [0] 


where 


’ ° = '3IV(to>’eco<to>'*c°<to)> 


nt,V  “ l3[b  ,b2,b3]  J 

'Jy'l 

. , s,  9 [right  hand  side  of  Eqns.  (17)] 
[ ( )]  = a[<t>co(t),0co(t)’,|'co(t)] 


(20a) 

(20b) 

(21a) 

(21b) 

(22a) 


3x3 


. , . , _ 9 [right  hand  side  of  Eqns.  (17)] 
1 U,J  ” 3[b1,b2,b3] 


(22b) 

The  partials  (22)  are  taken  by  direct  differentiation  of  Equations  (17). 


11 


2.4  OBSERVATION  EQUATIONS  FOR  DYNAMICALLY  CONSTRAINED  PHOTOGRAMMETRIC 
ADJUSTMENTS 


2.4.1  Basic  Notations 

Linearization  of  the  colinearity  equations  (1)  leads  to  the  matrix 
equation* 


AY  = A AX  + E (23) 

where  AY  = mxl  matrix  of  "measured  minus  computed"  residuals  of  x , 

Pij 


[A]  = [— J = m x n matrix  (m  > n)  of  partial  derivatives  of  the  (x  , 

Pij 

y ) with  respect  to  the  elements  of  the  vector  {X}  of 

Pij 

unknown  parameters,  evaluated  with  the  currect  best  estimates 
of  X.  (25) 

{AX}  = an  n x 1 matrix  of  desired  corrections  to  the  parameter 
vector  {X} 

{E}  = AY  - A AX  = an  m x 1 matrix  of  "errors  after  the  solution" 
associated  with  any  choice  for  AX.  (26) 


Subscript  ij  denotes  ith  point  measured  on  the  jth  photograph. 
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The  least  squares  solution  for  Ax  minimizes  EWE  and  is  well  known  to  be 

T -IT 

the  normal  equations  AX  = (A  WA)  A WAY  where  W is  a positive  definite 
weight  matrix. 

For  traditional  (unconstrained)  photogrammetric  adjustments,  the 
{X}  vector  can  be  partitioned  as 

!{XYZC}| 

He*}  [ (27) 

{XYZP}| 

where 


{XYZC}  = 


{ ipdip }“ 


<f>l 

0i 


= matrix  of  camera  exposure 
stations'  coordinates 


(28) 


matrix  of  "1-2-3"  camera  orientation  angles 
(depending  upon  the  option  used,  these  may  (29) 
relate  camera  axes  to  the  Instantaneous  radial 
transverse,  and  normal  directions;  or  directly 
to  the  earth-fixed  axes). 
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{XYZP  }* 


'l\ 


matrix  of  pass  points'  coordinates 


(30) 


and  the  matrix  of  partial  derivatives  can  be  similarly  partitioned  as 
[A]  = [ [AXYZC ] [A4>e*J  [AXYZP ] ] (31) 

where  the  submatrices  are 

• Y ; . . . ) 


[AXYZC] 


[A$e<j>] 


3(X  , Y , X , Y ; 

P11  P11  P21  P21 


Plj'  PiJ 


3(X  , Y > Z ; X , Y , Z >••••  X , Y , Z ;...) 

C-  C_  C-  '-o  '-"1  w 4 4 


L C1  C1  “1  ^2  *'2  v'2  ~J 


9(X  , Y ; X , Y X , Y ; ...) 


(32) 


Pll'  pll'  P21  P21 


9 (<(> -^ » '1*^5  4’2*  ®2'  ^2* * * * * ’ ' * * ^ 


(33) 


and 


[AXYZP ] 


3(X  , Y ; X^  , Y : 

P11  P11  P21  P21 


; X , Y ; ...) 

^ Plj_ 


8(Xpl’  Ypl‘  Zpl5  Xp2*  Yp2*  Zp2?  Xpi’  Ypi*  Zpi! 

If  one  introduces  translational  dynamical  constraints,  then  (for  a single 

strip  of  photographs  and  using  the  simplifying  assumption  that  the  n ex- 
posure times  are  perfectly  measurable)  the  {X}  vector  (27)  and  A matrix 
(31)  become 

{0E} 

(He*)  > (35) 


{X} 


{XYZP} 


(34) 
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and 


[A]  = [[AOE]  [A<f>0<H  [AXYZPJ] 


(36) 


where 


a set  of  six  orbit  elements  defining  the  best 

fitting  orbit  (for  example,  the  initial  orbit 

state  vector(X  , Y , Z , X , Y , Z ). 

o o o o o o 


[AOE] 


[3 (X  , Y ; X , Y ; 

P11  P11  P21  P21 
3(0E1,  0E2,  ....  0E6) 


x , y ; . • • ) 

pii  PU 


If  one  introduces  both  translational  and  rotational  dynamical  constraints, 
then  (again,  for  a single  strip  of  photographs),  the  {X}  vector  (27,35) 
and  the  A matrix  (31,36)  become 


(X)  = < <AE}  ) 
UXYZP}] 


[A]  = [[AOE]  [AAE]  [AXYZP ]] 


where 


[AAE] 


a set  of  six  attitude  elements  defining  the 
best  fitting  attitude  solution  (for  example, 
the  initial  orientation  state  vector,  , 0o> 

V bl‘  b2 ’ b3) 


~3  (X  , Y ; X , Y ; Xn  , Y ;...)“! 

P11  P11  P21  P21  Pi1  Pij 

_3(AElf  AE2 AE6)  J 


The  above  symbology  may  prove  difficult  to  extrapolate  to  specific 
situations,  we  therefore  provide  a simple  prototype  adjustment  problem 


I 
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as  an  illustration.  It  will  also  be  evident  that  multiple  strip  (block) 
adjustment  present  no  particular  difficulties;  independent  dynamical 
constraints  are  introduced  for  each  strip  of  photography. 


f 1 


I 


i 


L 


2.4.2  An  Example  Adjustment  Problem  Setup 

Referring  to  Figure  3,  an  idealized  block  adjustment  of  2 strips 
containing  5 photographs  each  is  considered.  Six  complete  ground  control 
points  are  indicated  by  19  pass  points  are  located  at  each 

and  under  each  photo's  principal  point  ("+").  Thus,  in  this  idealized 
adjustment  problem,  each  of  the  ten  photographs  has  at  least  one  control 
point;  photos  2,  3,  4,  7,  8,  and  9 have  a total  of  9 measurable  points 
(108  measurements),  while  photos  1,  5,  6 and  10  have  6 measurable  points 
(48  measurements) . 

We  first  consider  the  classical  approach  in  which  no  dynamical  con- 
straints are  imposed.  The  total  number  of  equations  (one  for  each  measur- 
ment)  is 

# of  photos 

2 E [measurable  points  in  photography]  = 2 [ 6+9+9+9+6+6+9+9+9+6 ] 

j=1  = 156  (42) 


The  number  of  unknowns  associated  with  the  19  pass  points  is 

//  of  elements  in  {XYZP}  = 3 x 19  * 57  (43) 

The  number  of  unknown  camera  exposure  station  coordinates  is 

# of  elements  in  {XYZC}  = 3 x 10  = 30  (44) 

The  number  of  unknown  orientation  coordinates  is 

//  of  elements  in  { 4>0^ } * 3 x 10  = 30  (45) 

Therefore,  the  total  number  of  unknowns  (without  dynamical  constraints) 


is  57  + 30  + 30  = 117. 
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COLUMN 


Fig.  3 An  Example  Triangulation  Block  of  2 Strips  of  5 Photo- 
graphs Each  (6  complete  ground  control  points  are  de- 
noted by  0,  the  10  camera  exposure  points  are  indicated 
by  +,  19  pass  points  are  located  at  the  • and  + symbols, 
50%  overlap  is  assumed  along  each  strip) 
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The  structure  of  the  [A]  matrix  for  this  example  (without  constraints) 
is  graphically  illustrated  in  Figure  4.  The  left,  middle,  and  right  parti- 
tions of  Figure  2 correspond  identically  to  the  partitions  of  Equation  (35). 
The  shaded  areas  are  the  nonzero  entries. 

We  now  consider  the  changes  which  occur  in  the  structure  of  the  partial 
derivative  matrix  when  translation  (orbital)  constraints  are  incorporated. 
Since  two  orbital  arcs  are  present,  then  one  dynamical  (best  fitting)  solu- 
tion will  be  recovered  for  each  strip.  Thus,  we  "trade"  the  30  unknown 
elements  in  the  {XYZC}  array  for  2 x 6 = 12  unknown  orbit  elements  in  a 
12  x 1 {OE}  array,  thereby  reducing  the  total  number  of  unknowns  to  117- 
18  = 99.  The  [A]  matrix,  now  modified  to  reflect  the  translational  con- 
straints is  shown  in  Figure  5.  Comparison  with  Figure  2 reveals  that  only 
the  left-most  partition  of  elements  is  modified,  but  that  considerable 
(60%)  compression  of  this  partition  has  been  achieved.  Not  considering 
exposure  times  to  be  perfectly  measurable  would  result  in  10  additional 
adjustable  parameters  (109  unknowns)  and  a 27%  compression  of  the  first 
partition  of  [A] . 

We  finally  consider  the  structure  when  both  translational  and  rota- 
tional constraints  are  imposed.  We  exchange  the  30  unknown  orientation 
parameters  for  2 x 6 = 12  unkonwn  orientation  elements  (1  set  of  six  for 
each  of  the  2 strips).  The  [A]  matrix,  sketched  in  Figure  6 now  shows  a 
60%  compression  of  the  center  partition  (the  total  number  of  unknowns  is 
now  117  - 36  = 81).  This  60%  compression  of  the  center  partition  is 
achieved  even  if  one  must  also  recover  the  n exposure  times,  since  these 
times  would  have  already  been  included  in  incorporating  the  orbital 
constraints.  In  practice,  however,  the  quality  of  the  gyro  data  and 
auxiliary  attitude  measurements  will  dictate  whether  or  not  this  theo- 
retical truth  results  in  improved  or  degraded  precision.  The  simulations 


Partial  Derivative  Matrix  [A]  for  Observation 
Equations  for  the  2 Strip  Block  for  Fig.  3 
WITH  TRANSLATIONAL  DYNAMICAL  CONSTRAINTS 


Fig.  6 Partial  Derivative  Matrix  [A]  for 
Observation  Equations  for  the  2 
Strip  Block  of  Fig.  3 

WITH  TRANSLATIONAL  AND  ROTATIONAL  DYNAMICAL 
CONSTRAINTS 


make  a weak  case  for  increased  accuracy,  but  a strong  case  for  increased 


computational  efficiency. 


2.5  Sequential  (Kalman  Filter)  Approach  to  Dynamically  Constrained 
Triangulation 

2.5.1  The  Expended  Kalman  F UteA  AlgosUXhm 

Given  a system  described  by  nonlinear  equations  of  the  form 
X = F(t,X)  + Vx  (42) 

with  measurements  modeled  by 

Y = G(X)  + V (43) 

y 

where  V and  V are  Gaussian  noise  vectors  with  zero  mean;  it  is  assumed 
x y 

V and  V are  uncorrelated, 
x y 

If  a sequence  of  measurements  are  available 
{t^,  Y^,  ^~2’  ^2’  ' ’ ' 1 tk’  ^k’  * " 

Then  the  extended  Kalman  estimation  algorithm  (see  ref.  10)  is  defined 
by  the  recursions 

*k+l  (tk+l)  = ^Sc^k+l*  + K(tk+1)  [Yk+l  “ G(Xk(tk+l))^  (44a)* 

where 

X.  (t,  = forward  *nonlinear)  integration  of  Equation  (42)  from 

W 

K“k+1>  ‘ Pk(tk+1>  »T<W  + H(tk+l)  P<Vl>  “T<tk+l>rl  («»> 

k+1  k+1 

= Kalman  Gain  Mat/Ux 

p^W  = $(tk+r  V pk(tk)  $T(tk+r  ck>  + Q (44c) 

Pk+l^k+l^  = X " covariance  matrix  = [I  - K(ck+1)  H^tk+1^  Pk^k+l^  (44d) 


H = [— ] = partials  of  the  measured  variables  with  respect  to  the 
state  vector  (44e) 

<f(t,tk)  = nxn  state  transition  matrix,  the  solution  of  the  differential 
equation 


4>(t,t.)  = [-rd  $(t,t  );  4>(t,  ,t  ) = I 


(44  f ) 


ft  A 

Interpretation  of  subscripts:  X.(t^)  = "the  best  estimate  of  X at 
t,  based  upon  the  first  j data  subsets". 
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Q = process  noise  covariance  matrix  variance  of  x - vector  due 
to  the  presence  of  the  noise  vector  in  Equation  (42) . 

2.5.2  Particularization  of  the  Extended  Kalman  Estimation  Algorithm 
to  the  Constrained  Triangulation  Problem 

The  above  algorithm  can  be  applied  to  the  constrained  triangulation 

problem  by  making  the  following  definitions: 


X = state  vector  = 


{XYZP }) 


{ XY  ZP } ; 


{Foe>  = 


(F  } 

ae 


{F  } 
{Foe} 
{0ae} 


-yx/r^  + x-perturbation | 
-yy/r^  + y-perturbation ' 
-yz/r  + z-perturbation 


Equation  (17) 


(46a) 


(46b) 


(47) 


Y^  = the  kth  particular  set  of  measured  image  coordinates  ^ X2 ( 
of  images  appearing  in  the  kth  photograph  )y9/ 


H(t^)  = all  pairs  of  Equations  (1),  one  pair  corresponding  to  each 
image  measured  in  the  kth  photograph. 

, 3H,  = f3H 3H|  3H 3H 3H , 

l3Xj  3 (x,y,  z)  3 (x,y , z)  3(<t>,6,^)  3(bx,by,bz)  3(XYZP)'J 


rlH 0 

l3(x,y,z) 


3(<M,'J>) 


3H 

3 (XYZP) 


[ , - ] = partials  of  Equations  (1)  w.r.t.(X  ,Y  ,Z  ) 

o V. X,  1 , L)  C C C 


[ , , , v; — rr]  = partials  of  Equations  (1)  w.  r . t . (<f> , 0 ,ij/) 
0 (<t>  ,0 » V) 


rM j 

1 3 (XYZP) J 


= partials  of  equations  (1)  w. r . t . (X^ ,Y^ ,Z^) 


[-]  = 
l3Xj 


({F  } 

\ oe 

{F“} 

\ _ae 

M 0 } 

{_0E} 

{_AE} 

{XYZP} 


3{F  } 3{F  } 

oe  oe 

3 { OE } 3{AE} 

3{^ae}  3{  ^~ae} 

3{0E}  3{AE} 


6x6 

■3{^oe} . 
' 3{0E}  ' 


3x3  3x3 

r 0 I , 


(49a) 


(49b) 


[G]  = ^5 
r 


(X2-r2/3)  XY 


(Y2-r2/e) 


(r  -r  /3). 


+ oribtal  (49c) 

perturbations 


3{F  } 

— — p-pr-  = 0 [zero  unless  orbital  perturbations  which  depend  upon 
^ attitude  drag,  are  included  in  the  force  model] 


3{F  } 

ae 
3 { OE } 


= [0]  + perturbations 


L 


3{F  } 3(1,9, <i>,b,b,b) 

r -1  = r — _ — 1 — =_i  = 

1 9{AE}  l8(4>,9,.J/,bx,by,bz)J 


d(i>,d,h  9(^,6, 4-) 

d(<p,6,<(i)  9(bx,by>bz) 

- n n - 


(49d) 


3 (<t> , 6 ,^) 


] = partials  of  Equations  (17),  w.r.t.  (<t>,0,^) 


[ 9 (b  r^b^~^b~-)]  = Partia^s  Equations  (17),  w.r.t.,  (b^.b  ,bz) 
' x*  y ’ z'  xyz 


Process  Noise  Matrix  (50a) 


Q0  =0 


°gyro 
QXYZP  = ° 


1 0 
0 1 


(50b) 


The  above  algorithm  is  applied  to  a dynamically  constrained  triangulation 
problem  in  section  4.0. 

It  should  be  noted  that  the  {XYZP}  vector  includes  only  the  particular 
coordinates  of  the  points  measured  on  the  k.  + 1^  photograph;  thus  the 
computer  program  requires  careful  "bookkeeping"  which  dynamically  alters 
this  array  and  the  associated  matrices  in  the  Kalman  filter  algorithm  to 
remain  consistent  with  the  (k  + l)*"'1  set  of  elements  of  {XYZP}. 

The  Kalman  algorithm  convex gence  history  (ref.  10)  leads  to 
the  conclusion  that  the  above  algorithm  should  be  interated  through 
all  measurements  at  least  two  or  three  times  to  obtain  the  best  estimates. 
In  the  absence  of  process  noise  (from  gyro  measurements),  Junkins  (10) 
has  demonstrated  numerical  equivalence  between  the  iterated  Kalman 
algorithm  and  the  batch  least  squares  algorithm.  The  presence  of  process 
noise,  however,  precludes  a rigorous  application  of  the  batch  least  squares 
algorithm  and  thereby  leaves  the  sequential  Kalman  algorithm  as  the 
only  practical  alternative. 
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3.0  BATCH  TRIANGULATION  EXAMPLES 


> 


3.1  Example  1:  Torque-Free  Case 

In  this  example,  we  consider  an  uncontrolled  satellite  moving  in 
a 100  mile  circular  orbit.  Two-body  orbital  dynamics  are  assumed.  The 
rotational  dynamics  is  assumed  governed  by  Equations  (12)  and  (13)  with 
L^  = Ly  = Lz  = 0 (i.e.,  a rigid  spacecraft,  subject  to  zero  external 
torque).  Equations  (12)  and  (13)  have  analytical  solutions  for  this  case, 
and  analytical  solutions  for  all  required  partial  derivatives  have  been 
derived  (see  Ref.  7 and  8).  For  purposes  of  illustration,  a short 
strip  of  photography  was  simulated;  corresponding  to  four  photographs 
made  by  a camera  fixed  in  a rigid  triaxial  satellite  from  a 100  mile 
circular  orbit.  Two-body  translational  motion  and  torque-free  rotational 
motion  were  assumed.  A single  strip  of  four  photographs,  each  with 
approximately  50%  overlap  was  generated.  These  photographs  covered  a 
ground  strip  of  approximately  100  miles  by  250  miles. 

Four  control  points  and  eight  pass  points  were  chosen  on  the  ground 
strip  and  distributed  in  such  a fashion  as  to  give  a total  of  60  obser- 
vation equations.  The  image  coordinates  were  modified  to  simulate  1 mil 
measurement  precision.  It  was  further  assumed  that  the  camera  focal  length, 
f,  the  principal  point  offset  (x0»yQ)  anc*  the  camera  exposure  times  were 
perfectly  known.  Initial  guesses  of  the  pass  point  coordinates  and 
camera  center  coordinates  were  corrupted  from  their  true  values  by  an 
average  of  2 miles.  Initial  guesses  on  camera  axes  Euler  angles  were 
corrupted  by  an  average  of  1.5°. 

Three  separate  least  squares  differential  correction  programs  were 
written  to  illustrate  the  impact  of  dynamical  constraints: 

UCPHOTO  (unconstrained  photogramme try ) : all  parameters  appearing  in 
the  observation  equations  were  treated  as  independent.  The  dimensions 


26 


of  the  A matrix*  were  60  (observations)  x 48  (parameters) , where 
48  = 3 times  the  number  of  pass  points  plus  3 times  the  number  of  camera 
center  locations  plus  3 camera  axes  Euler  angles  at  each  exposure  time 
= 3(8)  + 3(4)  + 3(4). 

0CPH0T0  (orbit  constrained  photogrammetry) : pass  point  coordinates,  camera 
orientation  at  each  time  and  orbit  initial  conditions  were  considered 
independent.  The  A matrix  was  dimensioned  60  (observations)  x 42 
(parameters),  where  42  = 3 times  the  number  of  pass  points  plus  6 orbit 
initial  conditions  plus  3 Euler  angles  for  each  photograph  = 3(8)  + 6 + 
3(4). 

DCPHOTO  (fully  dynamically  constrained  photogrammetry):  pass  point 
coordinates,  orbit  initial  conditions  and  attitude  initial  conditions 
were  considered  independent.  The  dimensions  of  the  A matrix  were  60 
(observations)  x 36  (parameters)  where  36  = 3 times  the  number  of  pass 
points  plus  6 orbit  initial  conditions  plus  6 attitude  initial  conditions 
= 3(8)  +6+6. 

The  structure  of  the  A matrix  for  these  three  cases  is  shown  in 
Figures  7a,  7b,  and  7c,  respectively.  The  shaded  region  contains  the  non- 
zero elements.  The  essential  results  of  these  numerical  tests  were 
summarized  in  Table  1.  For  this  simple  problem  the  best  results  were 
obtained  using  full  dynamical  constraints.  (Even  though  the  dimension 
of  the  A matrix  was  reduced  from  60  x 48  to  60  x 36,  the  computer  run 
time  was  not  decreased) . This  is  because  the  matrix  inversion  savings 
were  offset  by  the  relatively  elaborate  dynamical  model.  The  matrix 
calculations,  however,  dominate  real  world  triangulation  problems. 
Consequently,  substantial  cost  reductions  may  be  achieved  in  practical 
applications . 
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Max  pass  pt  error  Mean  pass  Std.  dev.  CDC  6400 

Program  after  convergence  point  error  of  all  central  processor 

(true  - converged)  (ft)  errors  (ft)  time  (sec) 

(ft) 

UCPHOTO  205.08  12.24  73.58  5.82 

0CPH0T0  93.45  -6.36  37.47  5.66 

DCPHOTO  -57.88  -1.60  25.54  6.86 

Table  1 A comparison  of  the  effects  of  dynamical 

constraints  on  photogrammetric  triangulation. 


3.2  Example  2:  A 10  Photo  Strip  (Via  Batch  Least  Squares  Triangulation) 

This  example  is  similar  to  the  DC  Photo  example  above,  except  in 
two  respects:  (1)  instead  of  assuming  zero  torque  and  using  the  analytical 
solution  of  Ref.  7,  the  Euler  differential  equations  were  integrated 
with  the  true  torque  history  specified,  (2)  the  number  of  photos  in 
the  strip  was  increased  from  4 to  10  to  provide  more  realistic  dimensions 
for  the  demonstration. 

138  simulated  measurements  were  calculated  for  a total  of  69 
photographic  images  of  10  control  points  and  25  pass  points.  The  simu- 
lated measurement  errors  were  initially  set  to  zero  to  verify  that  the 
converged  solutions  (from  a wide  variety  of  starting  estimates)  re- 
covered true  angles  and  position  coordinates  to  acceptable  precision  (all 
angles  were  recovered  to  within  1 arc  sec  and  all  camera  and  pass  point 
coordinates  recovered  to  within  .15  ft.).  In  particular,  Tables  2 and 
3 summarize  the  pass  point  coordinates'  starting  estimates,  converged 
values,  and  true  minus  converged  residuals  in  miles;  analogous  results 
are  recorded  in  Table  3 for  the  orbit  and  attitude  state  variables. 

This  example  required  approximately  one  minute  of  central  processor  run 
time,  over  half  of  which  was  spent  manipulating  and  inverting  matrices 
(as  opposed  to  integrating  the  equations  of  motion) . 


29 


This  example,  albeit  for  an  idealized  case,  demonstrates  clearly  the 
validity  of  the  dynamically  constrained  approach  for  the  case  that  the 
actual  torque  history  can  be  accurately  modeled.  However,  this  restriction, 
along  with  the  rigid  satellite  assumption  are,  in  fact,  major  obstacles 
in  practical  applications.  To  circumvent  the  necessity  of  modeling 
torques  and  making  idealizing  assumptions  regarding  satellite  rigidity, 
we  incorporate  angular  rate  gyro  measurements.  The  availability  of 
angular  rate  measurements  permits  the  luxury  of  not  modeling  torques 
or  body  flexibility  [in  order  to  correctly  generalize  Equations  (12)], 
rather,  one  directly  integrates  Equation  (17).  This  approach  leads  to 
two  algorithms  (batch  and  sequential) , which  are  applied  to  example 
cases  below. 
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TABLE  2 

Batch  Least  Squares 

Pass  Point  Solution  Summary  for  Example  2 


Point 

Number 


1 


Starting  Converged  Converged  Minus 

Value  (miles)  Value  (miles)  True  Residual  (ft) 


x 1.000000  0.000004 
y 10.500000  10.000006 
z 3962.400000  3963.399997 


+0.02 

+0.03 

-0.02 


2 


x -37.000000  -38.999998 
y 40.000000  39.000004 
z 3962.000000  3963.000002 


+0.01 

+0.02 

+0.01 


3 


x -12.000000 
y 42.000000 
z 3962.000000 


-10.000004 

41.999995 

3962.500025 


-0.02 

-0.03 

+0.13 


IJ 


25 


x 44.000000  45.000004 
y 431.000000  429.999995 
z 3941.000000  3940.49981 


+ .02 

- .03 

- .10 


Notes  on  error  statistics: 

•Mean  error  (over  all  25  pass  points)  - .0001  ft. 

•RMS  error  (over  all  25  pass  points)  - 0.05  ft. 

•Convergence  was  achieved  from  the  starting  values  shown  in  4 
differential  corrections. 

•For  brevity,  only  4 of  the  twenty-five  pass  point  coordinates  are 
shown. 

•Since  perfect  measured  values  were  used  in  this  example,  the  above 
small  residuals  simply  reflect  the  fact  that  (a)  the  algorithm 
has  been  correctly  inplemented,  and  (b)  adequate  precision  and 
convergence  tolerances  were  employed. 
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TABLE  3 

Batch  Solution  Summary  for  Initial  Orbit  and  Attitude 
State  Variables  (Example  2) 


Starting 


Converged 


Converged  Minus 


Value 

Value 

True 

Residual 

X 

0.004200 

miles 

0.003801 

miles 

.005 

ft. 

o 

yn 

0.002700 

miles 

0.000702 

miles 

.01 

ft. 

o 

z 

4062.181900 

miles 

4063.181898 

miles 

-.01 

ft. 

•o 

X 

-0.000160 

miles/sec 

-0.000060 

miles/sec 

.00 

ft/sec. 

•o 

¥n 

4.911200 

miles/sec 

4.851200 

miles/sec 

.00 

ft/sec. 

o 

z 

o 

-0.000100 

miles/sec 

-0.000160 

miles/sec 

.00 

ft/sec. 

♦o 

-0.000042 

rad. 

-0.000020 

rad. 

+1.0 

arc  sec 

0° 

0.000100 

rad. 

0.000000 

rad. 

0.0 

arc  sec 

0.000028 

rad. 

0.000030 

rad. 

-0.9 

arc  sec 

* 

-0.002100 

rad/sec 

-0.001100 

rad/sec 

0.00 

rad/sec 

6° 

0.000150 

rad/sec 

0.000100 

rad/sec 

0.00 

rad/sec 

?o 

Notes : 

0.000910 

rad/sec 

0.000990 

rad/sec 

-0.00 

rad/sec 

These  values  correspond  to  the  same  converged  solution  (after  4 
corrections)  as  the  pass  point  summary  of  Table  2. 


3.3  Example  3:  A Batch  Estimation  Example  Using  the 
V-OiecX  R cute.  Inte-gnjition  Approach 

This  example  is  very  similar  to  example  2,  except  the  rotational 

dynamical  model  consists  of  equation  (17)  [instead  of  equations  (12) 


and  (13)].  The  true  angular  rate  history  was  assumed  perfectly  measured 

by  rate  gyros  and  then  integrated  by  4-cycle  Runge-Kutta.  Instead  of 

the  Example  1 and  2 "attitude  elements"  (d>  (t  ),0  (t  ) fip  (t  ); 

co  o co  o co  o 

<Dx(to)  ,<Dy(to)  ,wz(to)  ) i we  successively  improve  instead  the  initial 
orientation  angles  { 4>  (t  ),  0 (t  ) ,<Ji  (t  )}  and  the  gyro  bias  parameters 
(b^by.b^).  The  numerical  results  were  identical  to  Example  2,  this 
test  serves  only  to  indicate  the  valididty  of  the  basic  idea.  This 
example  was  repeated,  using  perfect  gyro  data  but  imperfect  image 
coordinates  (with  1 mil  measurement  error  variance) . The  converged 
pass  points  differed  from  their  true  values  by  an  average  of  1.8  ft. 
with  a standard  RMS  error  of  28.9  ft. 

In  real  applications  the  effects  of  gyro  bias  and  noise  must  also 
be  considered.  Thus,  the  final  example  of  section  4.0  uses  simulated 
gyro  data  (including  bias  and  noise)  a sequential  (Kalman  filter) 
estimation  algorithm. 


4.0  EXAMPLE  4 SEQUENTIAL  TRIANGULATION  SUBJECT  TO  DYNAMICAL  CONSTRAINTS 

Based  upon  the  formulation  of  section  2.5.2,  a Kalman  sequential 
triangulation  algorithm  has  been  developed.  The  10  photo  strip  of 
examples  2 and  3 is  selected  as  a basis  for  discussion  and  comparison  of 
results.  This  particular  computer  program  involves  a significant  amount 
of  book-keeping  to  dynamically  restructure  the  various  matrices  to 
correctly  reflect  the  set  of  pass  point  coordinates  being  estimated  on 
each  cycle.  As  of  this  writing,  this  algorithm  is  still  under  develop- 
ment and  test;  it  is  anticipated  that  this  work  will  culminate  with 
successful  numerical  demonstrations  in  a few  months.  Although  this 
particular  portion  of  the  work  was  in  addition  to  the  central  objectives 
of  this  research  project,  the  authors  feel  carrying  it  to  completion 
will  be  a valuable  contribution.  This  work  will  thus  be  discussed  in  a 
subsequent  report,  upon  completion  of  the  computational  aspects. 
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5.0  CONCLUSIONS  AND  RECOMMENDATIONS 

This  report  details  several  approaches  to  the  dynamically  con- 
strained triangulation  problem.  The  formulations  are  applied  to  obtain 
numerical  results  for  several  cases.  These  results  support  the  following 
conclusions: 

(1)  The  most  attractive  method  studied  for  batch  triangulation  is 
the  direct  rate  integration  approaches,  as  discussed  in 
sections  2.2  and  3.3. 

(2)  The  sequential  (extended  Kalman  filter)  approach,  discussed 
in  sections  2.5  and  4.0,  appears  more  attractive  if 

(a)  pA0C&44  n04J>&  (due  to,  for  example,  integration  of 
noisy  gyro  measurements)  are  degrading  the  solution 
to  a considerable  degree,  and/or 

(b)  extremely  long  strips  of  photography  are  being  tri- 
angulated. 

We  recommend  that  future  refinement  development  of  triangulation 
software  include  dynamical  constraints  in  one  of  the  two  above  fashions. 
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APPENDIX  A 


COORDINATE  FRAME  NOMENCLATURE  AND  KINEMATICAL  RELATIONSHIPS 

Important  to  any  complex  angular  motion  problem  is  the  choice  of 
orientation  parameters  and  a clear,  systematic  nomenclature  convention. 

The  choice  of  orientation  parameters  is  critical  in  that  numerical 
difficulties  may  be  avoided  entirely  or  at  least  relegated  to  attitudes 
which  will  never  occur  in  the  problem  at  hand;  a self-evident  nomen- 
clature is  useful  (particularly  when  several  coordinate  systems  are 
defined)  not  only  to  readers  but  to  the  investigators  themselves. 

Throughout  this  research  effort,  Euler  angles  have  been  utilized 
for  coordinate  frame  orientation.  It  is  a well-known  fact  that  each  of 
the  twelve  possible  Euler  angle  sets  possesses  a kinematical  singu- 
larity. In  order  to  render  the  singularity  problem  as  harmless  as 
possible,  the  Euler  angle  sequence  for  a given  coordinate  frame  relative 
orientation  has  been  selected  so  that  encountering  the  singular  orien- 
tation is  highly  unlikely.  Future  efforts  could  profitably  involve  the  use 
of  a four  parameter  orientation  description  (such  as  Euler  parameters) 
so  as  to  avoid  the  singularity  problem  entirely. 

Unless  there  are  reasons  to  the  contrary,  the  1-2-3  Euler  angle 
sequence  will  be  employed  for  simplicity  of  discussion.  For  future 
reference,  the  3x3  transformation  matrix  [DE]  orienting  a frame  D 
characterized  by  the  orthonormal  vector  set 

{d}  = jd2 

(—3 

relative  to  a frame  E characterized  by 


(A.  1) 
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Cl^de 

S*de 

-S*de 

c*de 
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0 

de 

0 


0 -s< 

1 


de 

0 


^de-1 


0 

0 

C*de 

S*de 

“S*de 

C<^de 

where  ^g*®^’  and  are  the  1-2-3  Euler  angles  shown  in  Fig.  4 and 


c0^e  = cos®de>  s0^e  = sin0^e,  etc.  The  relationship  between  the  two 


frames  is  compactly  written  as  the  matrix  product 

W = [DE(<t>de,0de,*de)He}.  (A. 4) 

Obtainable  from  an  elementary  kinematic  analysis  of  the  frames  displayed 
in  Fig.  A4  is  the  following  matrix  expression  for  the  components  (along 
the  {d}  unit  vectors)  of  the  angular  velocity  of  frame  D with  respect  to 
frame  E: 


D/E 


D/E 


”cedec4,de 

S*de 

o“ 

-C0deS*de 

C*^de 

0 

_S0de 

0 

1 

Kde 


de 


(A. 5) 


-I  ' de 


Finally,  both  members  of  Eq.  (A. 5)  may  be  premultiplied  by  the  inverse  of 
the  3x3  matrix  on  the  right  hand  side  to  obtain  expressions  for  the  time 
rates  of  change  of  the  Euler  angles: 


Kde 


de  i 


C0 


de 


de 


"C*de 

~8*de 

0 

i d/e, 
“i 

1 D/E  i 

| l“l 

I 

C0deS^de 

cQdeC*de 

0 

) D/E 

j 2 

, • t“<»dc.*de)l  "2D/E 

F 

, (A. 6) 

-s0dec*de 

S0deC*de 

C0de_ 

l D/E' 
\ 3 

\ I D/E  1 

I \“3  1 

1 

Contained  in  the  above  discussion  are  sub-and  superscript  conventions 
which,  while  appearing  unnecessary  in  this  simple  context,  will  be  utilized  to 
great  advantage  in  conjunction  with  the  several  coordinate  frames  about  to 
be  introduced. 
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We  now  turn  to  the  first  of  the  specific  coordinate  frames  selected 
for  use  in  the  present  problem.  The  inertial  frame  N (see  Fig.  A. 2)  is 


-j 


non-rotating  and  geocentric;  it  forms  the  basis  for  the  on-board  star 
catalog  (if  on-board  star  sensing  capability  exists)  and  is  character- 
ized by  the  unit  vector  set 


n = 


-II 


1—2/ 


n_ 
i —3  i 


(A. 7) 


An  "orbit  inertial"  frame  O'  (see  Fig.  A. 2)  is  defined  by  the  unit 
vectors 

r 

-t  _ ~o 

2-3  " ~ 

o 

r x r 
, — o — o 

-2  ~ |r  x f I 


(A. 8) 


2-i  = °2  X -3’ 

where  the  subscript  o denotes  initial  conditions.  The  {o’}  unit  vectors 
are  projected  onto  the  {n}  unit  vectors  by  the  direction  cosine  matrix 
[O'N]  as 

{o’}  = [0'N]{n}  (A. 9) 


The  nine  constant  elements  of  [O'N]  are  calculated  from  the  orbit  state 
vector  by 


O'N 


31 


O'N 


32 


O'N 


33 


x 

o 

r 

o 


r 

o 


z 

o 

r 

o 


°’N2i  - r 

^2 

°'N22  - IT 

h 

°'n23  ■ r 


(A.1U) 
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°’N11 

" 0’N22°,N33  - 

0'N23°,N32 

°'N12 

" °,N23°'N31  ' 

0’N21°,N33 

°’N13 

- °,N21°,N32  - 

0’N22°,N31 

1 hl  " 

• • 

(y  z - z y ) 
o o o o 

h2  " 

(z  X - x z ) 
o o o o 

h3  " 

(x  y - y X ) 
o o o o 

h = 

2 2 

(hx  + u2  + h 

2.1/2 

3 ’ 

(A. 11) 


The  orbit  frame  0 (see  Fig.  A. 3)  rotates  about  o^*  c^e  orbit  normal, 

with  the  constant  angular  velocity  (relative  to  N or  O' ) 

0/N  0/0' 


“ n— 2’ 


where  ft  is  given  by 
2 ii 


ft 


2n  3/2 
a . 


(A. 12) 


(A. 13) 


(orbit  period)  ^ 

GM  is  the  earth  gravitational-mass  constant  and  a is  the  semi-major  axis 
of  the  orbit.  The  {<j}  unit  vectors  are  written  in  terms  of  the  {o'}  unit 
vectors  by 

(o)  = [00'(t)]{o'},  (A. 14) 

where 


[00* (t) ] 


cos[ft(t  - tQ)]  0 -sin[ft(t  - tQ)] 


sin[ft(t  - t )]  0 cos [ft(t  - t )]_J 

- o o 


(A. 15) 


Substitution  of  Eq.  (2.9)  into  Eq.  (2.14)  yields 

{o}  = [00' (t) ] [0'N] {n}  (A. 16) 

The  next  coordinate  frame  (G)  considered  is  one  attached  to  the 
strapped-down  gyroscopic  assembly  and  hence  body-fixed  but  not  nec- 
essarily principal.  The  G frame  is  generally  considered  to  be  the  pri- 
mary spacecraft  axes.  The  unit  vector  set  {g}  of  Fig.  A. 3 is  charac- 
terized by  being  the  roll  axis  (nominally  along  the  velocity  vector). 
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£2  being  the  pitch  axis  (nominally  along  the  orbit  normal) , and  being 
the  yaw  axis  (nominally  along  the  radius  vector).  The  1-2-3  angles  <t>gQ, 

0go>  and  i^o  orienting  G relative  to  0 are  hence  normally  small,  thereby 
avoiding  the  kinematical  singularity  at  0gQ  * 90°.  The  {g}  unit  vectors 
are  projected  onto  the  {o}  unit  vectors  by 

{&}  = [G0(4>go,ego.<l'go)]{°}-  (A- 17) 

Utilizing  Eq.  (A. 16): 

(g)  = [GO] [00 ’ ] [O'N] {n} . (A. 18) 

We  now  introduce  the  frame  C associated  with  the  terrain  camera  fields- 

of-view.  Frame  C is  characterized  by  c ^ being  along  the  f ield-of-view 

boresight,  with  and  £„  being  image  plane  axes.  The  angles  <J>  , 0 , 

i z eg  eg 

and  ip  orient  the  {c}  unit  vectors  relative  to  {g}  via  the  3-1-3  rotation 
eg  — a- 

{c}  = [CG(^cg,0cg^cg)]{&}  (A.  19) 

= [CG] [GO] {0}  (A. 20) 

= [CG][GO][00'J{0'}  (A. 21) 

= [CG] [GO] [00r ] [O'N] [n } . (A. 22) 

A 3-1-3  sequence  has  been  chosen  for  this  orientation  interlock,  with  the 
nominal  values 

^cg’W  = (0o,0°,0°).  (A. 23) 

The  orientations  prescribed  in  Eqs.  (A. 23)  will  generally  be  subject  to 
manufacturing  misalignments  and,  in  orbit,  will  likely  be  subject  to 
thermal  cycling  variations  of  a few  arc  minutes  (these  effects,  while 
important  in  applications,  are  not  treated  here;  we  note  that  their 
effects  can  be  absorbed  into  slightly  aliased  biases,  which  are  sequen- 
tially updated). 

Finally,  we  introduce  an  "orbit  C frame",  designated  C' , which  is 

coincident  with  frame  A at  the  initial  time  t , but  which  thereafter 

o 

rotates  about  the  orbit  normal  £2  * £2  with  the  constant  angular  velocity 
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C'/N  C'/O'  o « n / * 9 a \ 

w =w  = ft£2  = ft£2*  CA.  24) 

Frame  C'  represents  the  nominal  position  of  camera  frame  C and  is  oriented 

relative  to  frame  0 via 


{£'}  = c'0(<j>c,o,ec,o,ijic,o)  ] {£> , (A. 25) 

Where  the  1-2-3  Euler  angles  ((>  , , 0 , , and  4>  , are  constants  (nominally 

COCO  CO 

zero).  The  relationship  between  the  frames  C and  C'  is  given  by 

{c}  = [CC'(<f>ccl,eccl,*cc,  )]{£'},  (A. 26) 

wherein  the  three  1-2-3  angles  <|>  , , 9CC'»  and  ^cc*  are  normally  small. 

Table  A1  summarizes  the  coordinate  system  definitions  and  relationships 
discussed  thus  far. 

The  interlock  angles  (d>  ,0  , 4>  ) will  vary  about  their  nominal 

calibrated  values  due  to  thermal  cycling  or  any  other  structural  de- 
formation. The  maximum  amplitude  of  such  variations  is  believed  to  be  a 
few  arc-minutes.  We  have  not  yet  structured  the  estimation  algorithms 

to  provide  updates  of  the  interlock  angles  (<f>  ,0  , \p  ).  It  is 

eg  eg  eg 

believed  that  noise  in  typical  rate  data  (around  one  arc-sec/sec)  is  too 

large  to  allow  significant  refinement  of  (4>  ,0  , <p  ).  As  is  noted 

eg  eg  eg 

in  Section  4,  these  interlock  angles  are  correlated  with  the  estimated 
rate  bias  parameters,  a fact  which  also  hampers  accurate  estimation,  but 
also  makes  their  recovery  less  important;  small  interlock  angle  errors 
can  be  very  nearly  accounted  for  by  absorbing  these  errors  into  slightly 
aliased  rate  biases  (which,  in  themselves,  are  not  of  crucial  importance). 

Useful  for  consideration  in  incorporating  attitude  dynamical 
constraints  are  the  time  rates  of  change  of  various  Euler  angle  sets  in 
terms  of  the  inertial  angular  velocity  components  of  frame  G.  Speci- 
fically, we  will  develop  expressions  for 
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REFERENCE  FRAME  SUMMARY 


(1) 

(*  ,0  . 

) 

gn 

gn 

gn 

(2) 
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i. 

) 

cn 

cn 

cn 
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($  ,e  , 
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) 

go 

go 

go 

(4) 

(4>  t .0 
cc  cc 

.••*CC 

The  relationships  for  the  time  rates  of  change  of  the  angles 

f>  , 0 , and  ip  follow  directly  from  the  form  of  Eqs.  (A. 6): 

gnf  gn’  gn 

» “1/  G/N»  / G/N\ 

cili  -sty  0 liw,  1 |w, 

gn  gn 
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gn' 
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cB  si|>  c0  ci|/  0 
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„ s0  c^ 

*-  gn  gn 


s0  al>  c0  I \o) 
gn  gn  gn-1  ' 


L0G/N)  = [B(0  .*  ))  <*> 


G/N 


gn  gn 


G/N( 

2 

G/N  1 


(A 


Note  that  0 , for  the  case  of  a near  circular  orbit  and  one  vehicle 

gn 

rotation  about  Per  orbit»  will  go  to  90°  twice  per  orbit,  causing 

cos0gn  to  vanish  with  a resulting  singularity  in  Eqs.  (A. 27). 

The  relationships  for  the  time  rates  of  change  of  the  angles  <J>cn, 

0 , and  can  also  be  obtained  from  Eq.  (A. 6): 

cn  cn 


cn  y 

(0  \=  [B(0  ,4>  ,0  ,\f>  )]  <u)„ 

cn/  1 v an  an  ag’  ag,rag/J  ) 2 


♦ 


G/N 

L 

G/N 

G/N 


(A. 28) 


cn 


The  factor  of  [CG]  on  the  right  side  of  Eq.  (A. 28)  is  necessary  to  convert 

the  u>  from  components  along  {j^}  axes  to  components  along  {£>  axes.  Due 

to  the  form  of  the  matrix  B (defined  in  Eq.  (A. 6)  and  repeated  in  Eq. 

(A. 27),  singularities  would  occur  twice  per  orbit  for  a vehicle  which 

rotates  about  £2  once  per  orbit.  Hence,  although  fairly  simple  in  form, 

Eqs.  (A. 27)  and  (A. 28)  share  similar  difficulties. 

We  now  proceed  to  develop  equations  analogous  to  (A. 27)  and  (A. 28) 

for  the  angles  (<)>  ,0  , ip  ) (Case  3).  The  angular  velocity  of  G 

go  go  go 

relative  to  N may  be  written  as 


.27) 
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ooG/N  = uiC/O  + 10O/N 


(A. 29) 


• M}  • (2) 

a), G/N  + w„G/N  + oo.G/N  = <f>  o,  v ' + 0 o„v  ' + 4>  + fto0,  (A. 30) 

1 % 2 £2  3 % g0_1  gCr‘2  80  “2 
where  the  parenthetical  superscripts  indicate  intermediate  frames  in  the 

0 to  G transformation.  Expressing  o^3\  Oj^2^ , and  £2  *-n  terms  of 

^2»  and  enables  one  to  equate  coefficients  of  like  unit  vectors, 

resulting  in  three  scalar  equations  which  have  the  matrix  form 


0g°[  " [B(0go’V)1  p/N  | 

III  | G/N 

/ Iid, 

'go/  \ 3 


c<t>  si p + sd>  s0  cil i 

go  go  go  go  Tgo 


- 0 \C(t  ci li  - s4>  s0  s\l> 

) go  go  go  go  go 


-s<|>  C0 
go  go 


(A. 31) 


For  the  case  of  one  vehicle  rotation  about  Per  near  circular  orbit, 

the  angles  4>  , 0 , and  it  remain  small  so  that  Eqs.  (A.  31)  are  non- 
Ygo  * go’  ~ go 

singular  and  nearly  linear  (important  in  the  implementation  of  estimation 
algorithms) . 

• • • 

To  develop  the  Case  4 equations  for  <t>cc,,  ®cc'»  an<*  ^cc'»  we  begin 


G/N  C/C'  , C'/N 
0)  =0)  +0) 


(A. 32) 


Proceeding  in  a manner  analogous  to  that  outlined  for  Case  3,  we  arrive  at 
ll  \ .G/N, 


9 = [B(0  .,  * .)]  [CG(*  ,9  .*  ) /u) 


cc  cc 


eg  Cg  eg 


!cd>  si]/  + si  s0  ci 1/  \ 
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Like  Eqs.  (A. 31),  Eqs.  (A. 33)  are  more  complex  than  Eqs.  (A. 27)  or  (A. 28) 
but  affer  the  advantages  of  near  linearity  and  no  singularities  for  the 
case  of  one  vehicle  rotation  about  £2  per  orbit. 
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